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We present a study of different models of local disorder in graphene. Our focus is on the main 
effects that vacancies — random, compensated and uncompensated — , local impurities and sub- 
stitutional impurities bring into the electronic structure of graphene. By exploring these types of 
disorder and their connections, we show that they introduce dramatic changes in the low energy 
spectrum of graphene, viz. localized zero modes, strong resonances, gap and pseudogap behavior, 
and non-dispersive midgap zero modes. 
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Graphene is poised to become a new paradigm in solid 
state physics and materials science, owing to its truly bi- 
dimensional character and a host of rich and unexpected 
phenomena 1,2,3 . These have cascaded into the literature 
in the wake of the seminal experiments that presented 
a relatively easy route towards the isolation of graphene 
crystals 4 . 

Carbon is a very interesting element, on account of its 
chemical versatility: it can form more compounds than 
any other element 5 . Its valence orbitals are known to hy- 
bridize in many different forms like sp 1 , sp 2 , sp 3 , and 
others. As a consequence, carbon can exist in many 
stable allotropic forms, characterized by the different 
relative orientions of the carbon atoms. Carbon binds 
through covalence, and leads to the strongest chemical 
bonds found in nature. Common to the most interesting 
forms of carbon is the so-called graphene sheet, a single 
plane of sp 2 carbon organized in an honeycomb lattice 
(Fig. 1(a)). Graphite, for instance, is made of stack- 
ings of graphene planes, nanotubes from rolled graphene 
sheets, and fullerenes are wrapped graphene. Yet, for 
many years, it was believed that graphene itself would 
be thermodynamically unstable. This presumption has 
been overturned by a series of remarkable experiments 
in which truly bi-dimensional (one atom thick) sheets of 
graphene have been isolated and characterized 4 . This 
means that studies of the 2D (Dirac) electron gas can 
now be performed on a truly 2D crystal, as opposed to 
the traditional measurements made at interfaces as in 
MOSFET and other structures 6 . 

The crystalline simplicity of graphene — a plane of 
sp 2 hybridized carbon atoms arranged in a honeycomb 
lattice — is deceiving. The characteristics of the hon- 
eycomb lattice make graphene a half-filled system with 
a density of states (DOS) that vanishes linearly at the 
neutrality point, and an effective, low energy quasiparti- 
cle spectrum characterized by a dispersion which is linear 
in momentum 7 close to the Fermi energy. These two fea- 



tures underlie the unconventional electronic properties 
of this material, whose quasiparticles behave as Dirac 
massless chiral electrons 8 . Consequently, many phenom- 
ena of the realm of quantum electrodynamics (QED) 
find a practical realization in this solid state material. 
They include: the minimum conductivity when the car- 
rier density tends to zero 9 ; the new half-integer quantum 
Hall effect, measurable up to room temperature 9 ; Klein 
tunneling 10 ; strong overcritical positron-like resonances 
in the Coulomb scattering cross-section analogous to su- 
percritical nuclei in QED 11,12 ; the zitterbewegung in con- 
fined structures 13 ; anomalous Andreev reflections 14,15 ; 
negative refraction 16 in p-n junctions. 

Arguably, the most interesting and promising proper- 
ties from the technological point of view are its great crys- 
talline quality, high mobility and resilience to very high 
current densities 1 ; the ability to tune the carrier density 
through a gate voltage 4 ; the absence of backscattering 17 
and the fact that graphene exhibits both spin and valley 
degrees of freedom which might be harnessed in envis- 
aged spintronic 18,19 or valleytronic devices 20 . 

Disorder, ever present in graphene owing to its exposed 
surface and the substrates, is the central concern of this 
paper. In particular, we focus on the effects of vacan- 
cies and random impurities in the electronic structure 
of bulk graphene. The models examined below apply 
to situations in which Carbon atoms are extracted from 
the graphene plane (e.g. through irradiation 21 ), in which 
adatoms and/or adsorbed species attach to the graphene 
plane 22 , or in which some carbon atoms are chemically 
substituted for other elements. They are, therefore, mod- 
els of local disorder. We do not consider explicitly other 
sources of disorder like rough edges or ripples 23 , or the 
dramatic effects of Coulomb impurities, which have been 
discussed elsewhere 11,24 . In this article we expand the 
discussion of vacancies initiated in Ref. 25, using the 
same techniques, and discuss the consequences of local 
disorder originally presented in Ref. 26. Numerically we 
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FIG. 1: (Color online) In (a) a± and a<2 are the primitive 
vectors that define the WS unit cell highlighted as the dashed 
hexagon. The lattice parameter, a, is ~ 1.4 A. The first BZ of 
the associated reciprocal lattice is shown in (b) , together with 
the points of high symmetry T, M and the two nonequivalent 
K and K'. 



resort to exact diagonalization calculations and to the 
recursion method 27 ' 28 . The latter allows the calculation 
of the DOS and other spectral quantities for very large 
system sizes with disorder. In our case, the calculations 
below refer to honeycomb lattices with 4 x 10 6 carbon 
atoms, a size already of the order of magnitude of the 
real samples, if not larger for some experiments. 

The article is organized as follows. In section I we 
present the basic electronic properties of electrons in the 
honeycomb lattice, mostly to introduce the notation and 
the details relevant for the subsequent discussions. In 
section II we present our results regarding the different 
models of disorder. This section is subdivided according 
to the different models of disorder studied: vacancies in 
II A and II B, local impurities in II C and substitutional 
impurities in II D. The discussion of the results is kept 
within each subsection and the principal findings of this 
paper are highlighted in the conclusion, in section III. 



I. ELECTRONS IN A HONEYCOMB LATTICE 

Graphene consists of carbon atoms organized into a 
honeycomb lattice, bonded through covalence between 
two sp 2 orbitals of neighboring atoms (Fig. 1(a)). The 
graphene plane is defined by the plane of the sp 2 orbitals. 
The saturation of the resulting a bonding orbitals, leaves 
an extra electron at the remaining 2p z orbital per carbon 
atom. Ideal graphene has therefore a half-filled electronic 
ground state. 

The Bravais lattice that underlies the translation sym- 
metries of the honeycomb lattice is the triangular lat- 
tice, whose primitive vectors a\ and a2 are depicted in 



Fig. 1(a). One of the consequences is the existence of two 
atoms per unit cell, that define two sublattices (A and 
B in the figure): indeed, the honeycomb lattice can be 
thought as two interpenetrating triangular lattices. This 
bipartite nature of the crystal lattice, added to the half- 
filled band, imposes an important particle-hole symmetry 
as will be discussed later. 

The electronic structure of graphene can be captured 
within a tight-binding approach, in which the electrons 
are allowed to hop between immediate neighbors with 
hopping integral t ~ 2.7 eV, and also between next- 
nearest neighbors with an additional hopping t': 

H = -tJ2 4cj - a J2 c Uj + h - c - • (!) 

(iJ) <(i,j» 

The presence of the second term introduces an asym- 
metry between the valence and conduction bands, thus 
violating particle-hole symmetry. To emphasize the two 
sublattice structure of the honeycomb, we can write the 
Hamiltonian as 



H = - t ^ - t ^2 b i ai + S 

ieA,5 ieB,5 

-t' 4 a i+A - f 



(2) 



with operators and hi pertaining to sublattices A and 
B respectively. The vectors S connect atom i to its imme- 
diate neighbors, whereas the A connect atom i to its six 
second neighbors. Fourier transforming eq. (2) and in- 
troducing a spinor notation for the sublattice amplitudes 
leads to 



H 



e 2 {k) ei(fc) 
k Ui(k)* e 2 (k) 



, with 



(3) 

Since the spin degree of freedom does not play a role in 
our discussion other than through a degeneracy factor, 
it will been omitted, for simplicity. The functions ei(fc) 
and 62(h) read 



ei(fc) = -t£\ 



-iS.k 



-iA.k 



(4) 



£2(k) = — 2t f cos(V3k y a) 
'V3 



4t / cos ——k v a cos -k x a 



|ei(/c)| 2 =3t 2 + 2t 2 cos(V3k y a) 



(5) 



4t 2 cos I ^-k y a J cos (^k x a 



where 62(h) alone is the dispersion relation of a trian- 
gular lattice, and yield, after diagonalization of (3), the 
dispersion relations for graphene: 



E±(k)=e 2 (k)±\t\i 3 



e 2 (k) 



(6) 
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FIG. 2: (Color online) (a) Band structure along the symme- 
try directions of the reciprocal BZ of the honeycomb lattice, 
(b) Band structure of graphene (t f = 0) with the two bands 
touching at the K and K' points of the BZ. 
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The two bands E±(k) are represented in Fig. 2(b) in 
the domain k x , y G [— 7r, tt]. This unusual bandstructure 
makes graphene very peculiar with valence and conduc- 
tion bands touching at the Fermi energy, at a set of points 
at the edge of the first Brillouin zone, equivalent to the 
points K and K' by suitable reciprocal lattice transla- 
tions. Its low energy physics is dictated by the dispersion 
around those two inequivalent points, which turns out to 
be linear in k. In fact, expanding (5) around either 
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one gets the so-called K.p effective bandstructure 7 : 
E(K + q) = -3t' ±v F \q\ + 0(q 2 ) , 



ir.-« (7) 



(8) 



with a Fermi velocity, vp (hvp = vf = 3£a/2, and we 
take units in which h = 1 and t = 1). When t' = the 
dispersion is purely conical, as in a relativistic electron 
in 2D. For this reason, the two cones tipped at K and K' 
are known as Dirac cones. The low-energy, continuum 
limit of (2) is given by 



H = vp J d 2 r ip\r) a ■ pip(r) . 



(9) 



where ij)(r) is a two dimensional spinor obeying the Dirac 
equation in 2D 29 . 

Some quantitative aspects of graphene's band struc- 
ture (6) are plotted in Figs. 2 and 3. In panel 2(a) the 
band dispersion is plotted along the symmetry directions 
of the BZ indicated in Fig. 1(b), and in panel 3 the DOS 
for different values of the nearest-neighbor hopping, t' 
are plotted. Focusing on the particle-hole symmetric case 
(t = 0), it is clear that, besides the marked van Hove sin- 
gularities at E = ±t, the most important feature is the 
linear vanishing of the DOS at the Fermi level, a fact 
that is at the origin of many transport anomalies in this 
material 3,30 . 



FIG. 3: (Color online) DOS associated with eq. (6) for dif- 
ferent values of the next nearest neighbor hopping i! '. 



Particle-hole symmetry in this problem arises from the 
bipartite nature of the honeycomb lattice, and is a gen- 
eral property of systems whose underlying crystal lat- 
tice has this nature. When we have a bipartite lat- 
tice, the basis vectors of the Hilbert space can be or- 
dered so that, for any ket, \(p), the amplitudes in sub- 
lattice A come first. For example, if {0^, 4>\, • • • , 
are the Wannier functions for the orbitals in sublattice 
A, and {0^, . . . , 4>b} the ones in sublattice B, then 
our ordered basis could be {</>\, • • • , oV^ ; 0^, . . . , </>^}. If 
the Hamiltonian includes hopping only between near- 
est neighbors, this means that it only promotes itin- 
erancy between different sublattices. The stationary 
Schrodinger equation then reads, in matrix block form 
in the ordered basis, 



h A i 
Jab 



E 



(10) 



Expanding we get 



hAB ¥b 

h AB 



: Eif A 

E(f B 



(h\ B h AB ) <p B = E 2 ip B , (11) 



and therefore, if E is an eigenstate, so is — E. For a 
half-filled system, the elementary excitations around the 
Fermi sea can be thought, as usual, as particle-hole pairs. 
Since in that case Ef = 0, particles and holes have sym- 
metric dispersions. This is completely analogous to the 
situation found in simple semiconductors or semimet- 
als, although matters are slightly more complicated in 
graphene because there are two degenerate points, K and 
K' in the BZ. Thus there will be two families of parti- 
cle and hole excitations: one associated with the Dirac 
cone at K, and the other with the cone at K' \ like in a 
mult i- valley semiconductor. 



II. LOCAL DISORDER IN GRAPHENE 

Disorder is present in any real material, graphene being 
no exception. In fact, true long-range order in 2D implies 
a broken continuous symmetry (translation), which vio- 
lates the Hohenberg-Mermin- Wagner theorem 31,32 . So, 
by this reason alone, defects must be present in graphene 
and, in a sense, as paradoxical as it might sound, are pre- 
sumably at the basis of its thermodynamic stability. 

But the study of disorder effects on graphene is moti- 
vated by more extraordinary experimental results. One 
of them is the study undertaken by Esquinazi et al. 21 in 
which highly oriented pyrolytic graphite (HOPG) sam- 
ples were irradiated via high energy proton beams. As 
a result, the experiments revealed that the samples ac- 
quired a magnetic moment, displaying long range ferro- 
magnetic order up to temperatures much above 300 K. 
This triggered enormous interest, since the technologi- 
cal possibilities arising from organic magnets are many 
and varied. Furthermore, carbon, being the most cova- 
lent of the elements, has a strong tendency to saturate 
its shell in its allotropes, and is somehow the antithe- 
sis of magnetism. Besides the moment formation, it was 
found that the magnitude of the saturation moment reg- 
istered in hysteresis curves was progressively increased 
with successive irradiations. This is strong evidence that 
the defects induced by the proton beam are playing a 
major role in this magnetism. In this context the study 
of defects and disorder in graphene gains a significant 
pertinence. 

In the following paragraphs we will unveil some de- 
tails and peculiarities that emerge from different models 
of disorder applied to free electrons in the honeycomb 
lattice. 



A. Vacancies 

Vacancies are one of the defects more likely to be in- 
duced in the graphene structure by proton irradiation. 
A vacancy is simply the absence of an atom at a given 
site. When an atom is removed two scenarios are possi- 
ble: either the disrupted bonds remain as dangling bonds, 
or the structure undergoes a bond reconstruction in the 
vicinity of the vacancy, with several possible outcomes 33 . 
In either case, a slight local distortion of the lattice is 
expected. In the following discussion, however, it is as- 
sumed that, as first approximation, the creation of a va- 
cancy has the sole effect of removing the tt z orbital at 
a lattice point, together with its conduction band elec- 
tron. In this sense, the physics of the conduction band 
electrons is still described by the Hamiltonian (1), where 
now the hopping to the vacancy sites is forbidden. 



1. Vacancies and a theorem 

Vacancies have an interesting consequence when t' = 0. 
If the distribution of vacant sites is uneven between the 
two sublattices, zero energy modes will necessarily ap- 
pear. This follows from a theorem in linear algebra 34 
and can be seen as follows. Assume, very generally, that 
we have a bipartite lattice, with sublattices A and B (It 
can be any bipartite lattice like the square or honeycomb 
lattices in 2D, cubic in 3D, etc.), and that the number 
of orbitals/sites in A(B) is Na(Nb)- Just as we did be- 
fore, the basis vectors of the Hilbert space can always 
be ordered so that any ket, has the amplitudes on 
sublattice A appearing first, followed by the amplitudes 
on sublattice B: 

|*) = (^,^ S ) = (^,^,...,^;^,<A|,...,^ B ). 

(12) 

We now consider an Hamiltonian containing only nearest- 
neighbor hopping, plus some local energy (e^, e#) on each 
sublattice. The corresponding stationary Schrodinger 
equation will then be (in matrix block form that respects 
the ordering of the basis) 

17 17 \h\ B e B l N J \<p B ) \<Pb) 

(13) 

where 1m is the M x M identity matrix, Hab a Na x Nb 
matrix, and ifA (<Pb) a vector in a subspace of dimension 
Na (N b ). 

To analyze the spectrum we note that 

hcp B = (E -£ A )(fA , 14 x 
ft 1 " (pa = (E- e B )(PB 

which, from cross-substitution, implies that 

h^hifs = (E- e A )(E - e B )<PB • (15) 

If we call A 2 to the (non- negative) eigenvalues of ft 1 " ft, the 
spectrum of H is then 

E=^ ±] [(^) + V. (16) 

The symmetry about (sa + £ b)/2 simply reflects the 
particle-hole symmetry. 

2. Uncompensated lattices 

States of a peculiar nature should appear when the 
number of sites in each sublattice is different. With- 
out any loss of generality we take Na > Nb- Since the 
block Hab in (13) is a linear application from a vector 
space having dim (A) = Na-, onto a vector space B with 
dim(£?) = Nb, it follows from basic linear algebra that 

• rank(ft^ J e) = rank(ft^ B ) = Nb; 



• h<AB Pb — has no solutions other than the trivial 
one; 

• h\ B (fA = has non-trivial solutions that we call 

From the rank-nullity theorem, 

rank(Zi^) + nullity {h\ B ) = N A , (17) 

and hence the null space of h AB has dimension: 
nullity (h AB ) = Na — Nb- Consequently, there are states 
of the form 

l*°> = (^;0) 

in which (p\ satisfies h\ B yP A = 0, that are eigenstates of 
TL with eigenvalue ea- 




= (e A - e B )0 



(18) 



Furthermore, since nullity (h} AB ) = Na — Nb implies the 
existence of Na~Nb linearly independent ip A , this eigen- 
state has a degeneracy of Na — Nb- It should be stressed 
that a state of the form (<^;0) has only amplitude in 
the A sublattice. Therefore, we conclude that, whenever 
the two sublattices are not balanced with respect to their 
number of atoms, there will appear Na — Nb states with 
energy Ea, all linearly independent and localized only 
on the majority sublattice. In addition, one can modify 
sublattice B in any way (remove more sites, for instance) 
that these zero modes will remain totally undisturbed. 

We remark that in the above the details of the hopping 
matrix Kab were not specified and need not be. The re- 
sult holds in general, provided that the hopping induces 
transitions between different sublattice only, and that the 
diagonal energies are constant (diagonal disorder is ex- 
cluded). 



3. Zero modes 

The case with sa — £b — is of obvious relevance 
for us, since our model for pristine graphene does not 
include any local potentials. In this situation, the above 
results imply that introducing a vacancy in an otherwise 
perfect lattice, immediately creates a zero energy mode. 
Now this is important because those states are created 
precisely at the Fermi level, and have this peculiar topo- 
logical localization determining that they should live in 
just one of the lattices. 

Even more interestingly, it is possible to obtain the 
exact analytical wavefunction associated with the zero 
mode induced by a single vacancy in a honeycomb lat- 
tice. This was done by the authors and collaborators in 
Ref. 25, and will not be repeated here. We only mention 
that the wavefunction can be constructed by an appro- 
priate matching of the zero modes of two semi-infinite 
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FIG. 4: (Color online) Selected eigenstates in a graphene 
sheet with 80 2 atoms containing a single impurity at the cen- 
ter (black dot). Only the region near the vacancy is shown. 

(a) The eigenstate with energy closest, but different, to zero.. 

(b) The eigenstate with E — 0. (c) and (d) show the presence 
of two quasi- localized eigenstates even with t = 0.2£. 



and complementary ribbons of graphene, and that, in 
the continuum limit, the wavefunction of the zero mode 
introduced by one vacancy has the form 25 



iK'.r 



iK.r 



x + iy x — iy 



(19) 



The important point is that the amplitude of this state 
decays with the distance to the vacancy as ~ 1/r, and 
thus has a quasi- localized character, although strictly 
not-normalizable. And such quasi- localized state appears 
exactly at the Fermi level. 

Should another vacancy be introduced in the same sub- 
lattice, we already know that another zero mode will ap- 
pear. However, the nature of the two zero modes will 
depend whether the vacancies are close or distant. In 
the latter case, the hybridization between the two modes 
should be small on account of the 1/r decay, and we can 
expect two states of the form (19) about each vacancy 
site. Of course significant effects in the thermodynamic 
limit can only arise with a finite concentration of vacan- 
cies, and for such analysis we undertook the numerical 
calculations described next. 
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4- Numerical Results - Single Vacancy 

The first calculation is the numerical verification of the 
exact analytical result for the localized state in (19). For 
that, we consider the tight-binding Hamiltonian (1) and 
calculate numerically, via exact diagonalization, the full 
spectrum and eigenstates in the presence of a single va- 
cancy. For some typical results we turn our attention to 
Fig. 4. There we plot a real-space representation of some 
selected wavefunctions. This has been done by drawing 
a circle at each lattice site, whose radius is proportional 
to the wavefunction amplitude at that site, and whose 
color (red/blue) reflects the sign (+/-) of the amplitude 
at each site. Thus bigger circles mean higher amplitudes. 
In the first panel, 4(a), we are showing the eigenstate 
with lowest, yet non-zero, absolute energy. It is visible 
that the wavefunction associated with such state spreads 
uniformly across the totality of the system, like a plane 
wave. In the second panel, 4(b), we draw the wavefunc- 
tion of the state E = 0, that corresponds to (19). The 
state is clearly decaying as the distance to the central va- 
cancy increases. In addition, the state exhibits the full C3 
point symmetry about the vacant site, just as expected. 
This picture provides a snapshot of the lattice version 25 
of (19). Since only one vacancy was introduced, the state 
shown in Fig. 4(b) is the only zero mode present. 

When particle-hole symmetry is disturbed by a non- 
zero t', we still find states having this quasi- localized na- 
ture, where the wavefunction amplitude is still quite con- 
centrated about the vacancy. Two examples are shown 
in panels 4(b) and 4(c). They are two eigenstates with 
neighboring energy calculated for the same system. An 
important difference occurs here, in that, unlike the 
case t' = where only one localized state appears, the 
particle-hole asymmetric case opens the possibility for 
more than one of such states. 

This fact can be seen more transparently through the 
inverse participation ratio (IPR) of the eigenstates. With 
such purpose in mind, the IPR 



V(E n ) = J2\*n(n)\ 4 



was calculated across the band in both the t' = and 
t' ^ cases, with a single central vacancy. Typical re- 
sults are shown in Fig. 5. From 5(a) we do confirm that, 
when t' = 0, the presence of a vacancy introduces a lo- 
calized state at E = 0, which is reflected both by the 
enhanced IPR there, and by the sharply peaked LDOS 
calculated at the vicinity of the vacancy site. Although 
not shown in this figure, the amplitude of the peak in the 
LDOS at E = 0, Pi(0), decays as the distance between 
Ri and the vacancy increases, in total consistence with 
the analytical picture. When next-nearest neighbor hop- 
ping is included, we also confirm the appearance of states 
with a considerably enhanced IPR. Not only that, but, 
instead of one, we do observe a set of states with IPR 
much larger than the average for the remainder of the 
band. The LDOS is also enhanced near these energies, 
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FIG. 5: (Color online) IPR and LDOS calculated at one site 
closest to the vacancy. In panel (a), we have results for the 
IPR with i! = and t ^ without any vacancy (top row), 
and with a single vacancy (bottom) for comparison. In panel 
(b) we show the dependence of the IPR of the zero mode, 
V(E = 0), with the system size N (left), and also (V(E)) 
versus N for the remainder (extended) states (right). Dashed 
lines are guides for the eye. 



although the effect appears as a resonance on account of 
the finite DOS, in contrast with the sharp peak in the 
previous, particle-hole symmetric, case. 

A more definite and quantitative analysis is provided 
by the results in the subsequent panel (Fig. 5(b)). Here 
we present the dependence of V(E) on the number of car- 
bon atoms in the system, N. To understand the differ- 
ences, we recall that the IPR for extended states should 
scale as 



(20) 



But, for the zero mode (it should be obvious that when 
the term zero mode is employed, we are referring to the 
case with t' = 0.), we face an interesting circumstance. 
Remember that the wavefunction (19) is not normaliz- 
able. So, strictly speaking, the state is not localized, 
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and hence the designation quasi-localized that we have 
adopted above. The consequence of this is that the nor- 
malization constant for Sfr(x,y) depends on the system 
size: 



N 



V)\ 2 ~ log(ViV) ~ log(JV) . (21) 



This, in turn has an effect on the IPR because V(E) is 
defined in terms of normalized wavefunctions: 



N 



V(0) = 



log(JV): 



■£l*(*>y)l 4 



log(AT)2 



(22) 



This scaling of the IPR with N is precisely the one ob- 
tained numerically in Fig. 5(b) (left) for the zero mode, 
and is just another way of confirming the 1/r decay of 
this wavefunction. 



5. Numerical Results - Finite Concentration of Vacancies 

Unlike the single vacancy case, the dilution of the hon- 
eycomb via the introduction of a finite concentration of 
vacancies is not solvable using the analytical expedients 
employed in ref. 25, and numerical calculations become 
essential in this case. Our procedure consists in diluting 
the honeycomb lattice with a constant concentration of 
vacancies, which we call x (x = N vac /N). The diluted 
sites are chosen at random and the global DOS, aver- 
aged over several vacancy configurations, is calculated 
afterwards. This is clearly a disordered problem, and we 
employ the recursive method allowing us to obtain the 
DOS for systems with 2000 2 sites (which is already of 
the order of magnitude of the number of atoms in real 
mesoscopic samples of graphene studied experimentally) . 
Some results are summarized in Fig. 6. One of the ef- 
fects of this disorder is, as always, the softening of the 
van-Hove singularities (not shown). But the most sig- 
nificant changes occur in the vicinity of the Fermi level 
(Fig. 6(a)). In the presence of electron-hole symmetry 
it' = 0), the inclusion of vacancies brings an increase of 
spectral weight to the surroundings of the Dirac point, 
leading to a DOS whose behavior for E w mostly re- 
sembles the results obtained elsewhere within coherent 
potential approximation (CPA) 30 . Indeed, for higher di- 
lutions, there is a flattening of the DOS around the cen- 
ter of the band just as in CPA. The most important 
feature, however, is the emergence of a sharp peak at the 
Fermi level, superimposed upon the flat portion of the 
DOS (apart from the peak, the DOS flattens out in this 
neighborhood as x is increased past the 5 % shown here). 
The breaking of the particle-hole symmetry by a finite 
t' results in the broadening of the peak at the Fermi en- 
ergy, and the displacement of its position by an amount 
of the order of t' . All these effects take place close to the 
the Fermi energy. At higher energies, the only deviations 
from the DOS of a clean system are the softening of the 
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FIG. 6: (Color online) IPR and DOS for the diluted hon- 
eycomb lattice, (a) DOS for selected concentrations, x, and 
different values of t' . (b) IPR for selected values of t' using a 
concentration x = 0.5%. For comparison, the corresponding 
DOS is also plotted in each case. The concentration of vacan- 
cies is x, and only the vicinity of the Fermi level is shown. 
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van Hove singularities and the development of Lifshitz 
tails (not shown) at the band edge, both induced by the 
increasing disorder caused by the random dilution. The 
onset of this high energy regime, where the profile of the 
DOS is essentially unperturbed by the presence of va- 
cancies, is determined by e « vp/l, with I ~ n im p being 
essentially the average distance between impurities. 

To address the degree of localization for the states near 
the Fermi level, the IPR was calculated again, via exact 
diagonalization on smaller systems. Results for different 
values of t' are shown in Fig. 6(b) for random dilution at 
0.5 %. One observes, first, that V m ~ S/N for all energies 
but the Fermi level neighborhood, as expected for states 
extended up to the length scale of the system sizes used 
in the numerics. Secondly, the IPR becomes significant 
exactly in the same energy range where the DOS exhibits 
the vacancy-induced anomalies discussed above. Clearly, 
the farther the system is driven from the particle-hole 
symmetric case, the weaker the localization effect, as il- 
lustrated by the results obtained with t' = 0.2 1. To this 
respect, it is worth mentioning that the magnitude of the 
strongest peaks in V m at t' = and t' = 0.1 1 is equal to 
the magnitude of the IPR calculated above for a single 
impurity problem. Such behavior indicates the existence 
of quasi-localized states at the center of the resonance, 
induced by the presence of the vacancies. For higher dop- 
ing strengths, the enhancement of Vm is weaker in the 
regions where the DOS becomes flat. This explains the 
qualitative agreement between our results and the ones 
obtained within CPA in that region, since CPA does not 
account for localization effects. 

In summary, in this section, we saw that a single va- 
cancy introduces a quasi- localized zero mode. Its pres- 
ence is ensured by the uncompensation between the num- 
ber of orbitals in the two sublattices, and a theorem from 
linear algebra. The presence of this mode translates in 
the appearance of a peak in the LDOS near the vacancy, 
and in an enhanced IPR for this state. When we go from 
one to a macroscopic number of vacancies, we saw that 
both the peak and the enhancement of the IPR persist 
in the global DOS at Ep . 



B. Selective Dilution 

It is important to recall that the results of the previous 
section pertain to lattices that were randomly diluted. 
During such process, we expect the number of vacancies 
in sublattice A to be equal to the number of vacancies 
in sublattice B, on average. Strictly speaking, since our 
original lattices are always chosen with Na = Nb-, the 
fluctuations on the degree of uncompensation, Na — Nb, 
should scale as thus vanishing in the thermody- 

namic limit. Because of this, in principle, we would ex- 
pect the lattices used above to be reasonably compen- 
sated. But the theorem in § II A 1 only guarantees the 
presence of zero modes when the lattice is uncompen- 
sated. It turns out that, notwithstanding our utilization 




FIG. 7: (Color online) Dilution of just one sublattice of the 
honeycomb, (a) DOS for different dilution strengths, diluting 
only sublattice A. (b) On the left we show a detail of the 
DOS and the evolution of the gap with vacancy concentra- 
tions. On the right we plot the dependence of the missing 
spectral weight on the band (= 1 — w$) with x (circles). The 
continuous line is the best fit using f(x) — ax/(b — x) to the 
data represented by the circles. 

of rather large system sizes, such y/N fluctuations are 
still significant and the lattices were indeed slightly un- 
compensated. 

This clearly begs the clarification of the origin of the 
zero modes in the cases with finite densities of vacancies. 
Do they appear only through these fluctuations in the 
degree of sublattice compensation, or can we have zero 
modes even with full compensation? To try to elucidate 
this we developed a controlled approach to this issue, in 
the following. From now on we consider only the particle- 
hole symmetric situation (t f — 0). 



1. Complete uncompensation 

We have studied the DOS for systems in which only 
one of the sublattices was randomly diluted, with a finite 
concentration of vacancies. In this case, the system has 
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FIG. 8: (Color online) The gap estimated from the numerical 
curves in Fig. 7 is plotted against the vacancy concentration, 
x. The continuous line is a least squares fit to f(x) = ax b . 



precisely a number of zero modes that equals the number 
of vacancies. Starting from a clean lattice with N = Na + 
Nb sites, the latter corresponds to N v = Nx. We should 
thus expect a 5(E) peak contributing to the global DOS, 
with an associated spectral weight, w$ that coincides with 
the fraction of zero modes: 



w s (x) 



Nx 



N(l - x) 



1 



(23) 



Since the total spectral weight is normalized to 1, the 
spectral weight at E = has to be transferred from the 
states in the band. In Fig. 7 we show what is happen- 
ing. As seen in 7(a), the selective dilution promotes the 
appearance of a gap in the DOS, whose magnitude in- 
creases with the amount of dilution. At the center of the 
gap we can only see an enormous peak (not visible in the 
range used) staying precisely at E = 0, corroborating our 
expectations regarding the Dirac-delta in the DOS. But 
since it appears exactly at E = 0, we cannot resolve nu- 
merically its associated spectral weight. To obtain such 
spectral weight we calculated the spectral weight loss in 
the remainder of the band. The result and its varia- 
tion with the amount of dilution, x, is displayed in the 
right-most frame of Fig. 7(b). A non-linear fit to the 
data reveals that the dependence expected from (23) is 
indeed verified by the accord between the fitted curve in 
Fig. 7(b) and eq. (23). 

As Fig. 7 shows, the spectral weight is transferred al- 
most entirely from the low energy region near Ep and 
from the high energy regions at the band edges. This 
depletion near E = introduces the gap, 2E g . A gap 
implies the existence of a new energy scale in the prob- 
lem. Since the hopping t is the only energy scale in the 
Hamiltonian, such new scale has to come from the con- 
centration of vacancies. By dimensional analysis, such 
scale is dictated essentially by the average distance be- 
tween vacancies (I) 
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FIG. 9: (Color online) DOS for the honeycomb lattice using 
the controlled selective dilution discussed in the text, calcu- 
lated for different concentrations of vacancies, x, and several 
degrees of uncompensation, r\. Only the low energy region 
close to the Dirac point is shown. 



When the the magnitude of the gap found numerically is 
plotted against x we arrive at the curve of Fig. 8. The 
least squares fit shown superimposed onto the numeri- 
cal circles leaves confirms this assumption, and we arrive 
at a quite interesting situation, of having a half-filled, 
particle-hole symmetric and gapped system, with a finite 
concentration of (presumably quasi- localized) zero modes 
at the mid-gap point. 



2. Controlled uncompensation 

We now turn to a more controlled approach to the di- 
lution and uncompensation. For that we introduce an 
additional parameter, 77, that measures the degree of un- 
compensation. As before, we want to study finite con- 
centrations of vacancies. This is determined by x in such 
a way that the number of vacancies in a lattice with TV 
sites will be N v = Nx. But now, the number of vacancies 
in each sublattice is determined by 

1 



Nf = ^Nx(l + r,) 

n* = \nx(\ - n) , 



(25) 



vacancies 



(h = 1) . (24) 



with < 77 < 1. Therefore, the parameter 77 permits 
an interpolation between completely uncompensated di- 
lution (77 = 1), and totally compensated dilution (77 = 0). 
Let us look directly at the results for the DOS, calculated 
at different x and 77, and plotted in Fig. 9. 

At any concentration x the following sequence of events 
unfolds as 77 decreases from 1 to 0: (i) There is a perfectly 
defined gap in the limit 77 = 1.0 discussed above; (ii) 
for 77 < 1 a small hump develops at the same energy 
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scale of the previous gap; (iii) although the gap seems 
to disappear, it is clearly visible that when 77 < 1, the 
DOS decays to zero after the hump in a pseudo-gap like 
manner, and is zero at E = 0; (iv) decreasing further 77 
towards complete compensation (say, for 77 = 0.6, 0.4), 
this behavior persists, being visible that the DOS drops 
to zero at E = 0; (v) closer to full compensation (77 = 0.1) 
the DOS seems to display an upward inflection near E = 
0, and apparently does not drop to zero. Unfortunately, 
we are unable to resolve this region numerically with the 
desired accuracy. For instance, at higher dilutions x = 
0.2 we can still see the curve of 77 = 0.2 dropping to zero 
near E = 0. 

Naturally that, for all the cases with 77 7^ 0, the ex- 
istence of — Ny zero modes is guaranteed. As be- 
fore, we inspected this by calculating the missing spectral 
weight in the bands, and confirmed that it does agree 
with the fraction of uncompensated vacancies. Hence, 
the picture emerging from these results seems to suggest 
that, although the gap disappears for 77 < 1, the DOS still 
drops to zero at E = 0, and might drop in a singular way 
as 77 approaches zero. If we separate the contributions 
of the zero modes to the global DOS from the contribu- 
tion of the other states, the consequence of this would 
be that, in a compensated lattice (77 = 0), the DOS as- 
sociated with the other states would seem to diverge as 
E — > 0, but would be zero precisely at E = 0. Stated 
in another way, coming from high energies, we would see 
a decreasing DOS up to some typical energy e ~ y/x, at 
which point it would turn upwards. At very small ener- 
gies the DOS would seem to be diverging but, at some 
point arbitrarily close to E = 0, it would drop precipi- 
tously down to zero. Unfortunately, at the moment the 
numerical calculations are not so accurate as to allow the 
confirmation or dismissal of such possibility. In fact, the 
peaks for 77 = 0.0 are of the same magnitude of the ones 
found when the dilution is completely random across the 
two sublattices (Fig. 6(a)). So, although the evidence is 
compelling towards the affirmative, these results are still 
inconclusive as to whether the zero modes disappear in 
a perfectly compensated diluted lattice or not. 

C. Local Impurities 

Vacancies are local scatterers in the unitary limit. A 
vacancy can be thought as an extreme case of a local po- 
tential, U, when U — > 00. In this context we investigated 
the intermediate case characterized by a finite local po- 
tential. The Hamiltonian in this case changes from the 
pure tight binding in (1) to 

The first term represents the local potential of magnitude 
U at the impurity sites p. These impurity sites belong 
to the underlying honeycomb lattice but their space dis- 
tribution is random. The concentration of impurities, 
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FIG. 10: (Color online) Real part of (30) (top) obtained from 
the homogeneous DOS (bottom) through the Kramers- Kronig 
relations. 

x = Ni/N, is kept constant and we consider only the 
case with t' = in the sequel. 

Physically the model summarized in the Hamilto- 
nian of eq. (26) could describe the situation in which 
some of the carbon atoms are substituted by a different 
species. Another realistic circumstance has to do with 
the fact that a real graphene sheet is expected to have 
some molecules from the environment adsorbed onto its 
surface 22 . Consequently, even if the honeycomb lattice 
of the carbon atoms is not disrupted with foreign atoms, 
the presence of adsorbed particles can certainly induce 
a local potential at the sites where they couple to the 
carbon lattice. 

Much of the details of this model can be understood 
from the local environment around a single impurity, in 
which case exact results and closed formulas are obtain- 
able within a T-matrix approach 35 . Hence, we start by 
analyzing the single impurity problem in the honeycomb 
lattice, taking into account the full electronic dispersion 
and calculating the exact local Green's functions, which 
allow the identification of the main spectral changes in- 
troduced by the scattering potential. Within T-matrix, 
the 2x2 electron Green's function is written as 

G r y = G Q rr , + ^2 Gr,sT s ,s'G° s , r , . (27) 

s 7 s' 

In the Dyson-like expansion above G° is the non- 
interacting Green's function whose matrix elements are 
denoted by [G ]^f, (sub/superscripts refer to posi- 
tion/sublattice), and the T-matrix, T, is formally defined 
in terms of the scattering potential, V, by 35,36 

T{E) = T=^v ■ (28) 

Taking V r a f, = U5 r yS rt o^a,/3^a,o f° r a potential localized 
only on site r = of sublattice A, the local Green's func- 
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FIG. 11: (Color online) (Top) LDOS at the impurity site 
(po(E)) for different strengths of the scattering potential, in- 
dicated in the legend. (Bottom) The same data divided by 
the free DOS (p°(E)). 



tion on that site reads 



G\ 



AA 
00 



[G 



00 



U[GXo A 



(29) 



The function [G ]^ is simply related to the density of 



states per carbon atom in the absence of impurity, p°(E), 
through 



[G°(E)}^ A =F(E)-inp°(E). 



(30) 



The knowledge of p° (E) suffices for the determination of 
F(E) on account of the analytical properties of G°(E) 
and the Kramers-Kronig relations. Moreover, any new 
poles of the exact Green's function can come only from 
the denominator in (29), and are determined by the con- 
dition 



1 - UF(E) = . 



(31) 



Should this condition be satisfied for E within the branch 
cut of G°, the new poles will signal the existence of res- 
onant states in the band, and bound states of the local 
potential otherwise. Since p°(E) is known exactly 37 (cfr. 
Fig. 3), so can be G°(E) through eq. (30). The function 
F(E) is shown in Fig. 10. The profile of this function 
and the condition above, allows two immediate conclu- 
sions without further calculation: (i) the presence of the 
local potential induces bound states beyond the band 
continuum and (ii) a resonance appears at low energies 
beyond a certain threshold, C/ res , with energy of opposite 
sign with respect to the scattering potential, and which 
moves toward E = as U — > oo. 

The latter characteristic is certainly more interesting 
and we explore it a little further. To that extent notice 
that from (29) and (30) follows the interacting LDOS at 
the impurity site: 



Po(E) 



P°(E) 



[(l-UF(E)Y+[nU(P(E)] 



2 ' 



(32) 
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FIG. 12: (Color online) (a) Position of the maximum in the 
LDOS compared with the roots of (31). (b) Energy of the 
bound state. 



This quantity was calculated using the results of Fig. 10 
and (32) and, at the same time, using the recursive 
method that we have been using so far. The two do 
coincide, just as expected since the solution of the sin- 
gle impurity problem is exact, and, on the other hand, 
the recursive method is exact for the particular case of 
the LDOS 27 . The LDOS at the site of the impurity is 
shown in the top frame of Fig. 11 for several values of U. 
The bottom frame shows the same data divided by the 
non-interacting DOS, which amounts to replacing p°(E) 
by unity in the numerator of eq. (32). The resonance al- 
luded above is visible in both panels through the marked 
enhancement of the LDOS in the vicinity of the Dirac 
point. The position of the maximum in po(E) differs 
slightly from the roots of eq. (31) due to the modula- 
tion introduced by p°(E) in (32). This effect is shown 
in detail in Fig. 12(a) where the two values are explicitly 
compared. In addition, the LDOS also exhibits the Dirac- 
delta peak associated with the bound state (not shown 
in the figure), whose energy is plotted in Fig. 12(b) as 
a function of U . It is worth mentioning that analytical 
expressions can be obtained for the resonant condition 
(31) using the low energy Dirac approximation to the 
electronic dispersion 30 ' 38 . 

Returning now to our initial goal of populating the 
lattice with a finite concentration of local impurities, we 
expect the main features of the above analysis to hold to 
a large extent. But new features should also emerge from 
the possibility of multiple scattering and interference ef- 
fects in a multi-impurity environment. Although some 
of these effects can be captured within standard approx- 
imations to impurity problems 36 , we choose to present 
the exact numerical results obtained with the recursion 
technique. Examples of such calculations are shown in 
Fig. 13, where the global DOS averaged over several con- 
figurations of disorder is shown for different potential 
strengths and concentrations. 

The presence of the local term clearly destroys the 
particle-hole symmetry, leading to the asymmetric curves 
in the figure. As Fig. 13(a) makes clear, among the fea- 
tures seen locally for a single impurity (Fig. 11), the ones 
that carry to the global DOS of the thermodynamic sys- 
tem with a finite concentration of impurities are the res- 
onant enhancement of the DOS in the vicinity of the 
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FIG. 13: (Color online) DOS of the honeycomb lattice with 
a finite density of local impurities, (a) shows the DOS for 
U — 12t and different concentrations of impurities (notice 
the truncation in the horizontal axis), (b) shows a detail of 
the low energy region for different U and x as noted in the 
different graphs. 



FIG. 14: (Color online) (a) Variation of the "Dime point" 
energy Ed with impurity concentration and strength. The 
inset shows Ed/EU as a function of U, in which the curves 
with U < 3 roughly collapse onto each other, (b) Spectral 
weight transfer to the impurity band in the presence of local 
impurities. The values shown correspond to integration of the 
global DOS beyond E = 4t (above the main band edge, cfr. 
Fig. 13). 



Dirac point, and the high energy features that dominate 
beyond the band edge, and are associated with the im- 
purity states. One verifies that a finite concentration, x, 
generates a sort of impurity band at scales of the order of 
U, in accordance with the results in Fig. 12(b). This im- 
purity band has an interesting splitted structure as can 
be seen in the figure and is completely detached from the 
main band for U > 4t. In Fig. 13(b), we amplify the low 
energy region and display what happens as U and x vary. 
At small x and U the the DOS changes only through a 
simple translation of the band with the concomitant shift 
in the Dirac point, Ed- This rigid shift of the band at 
low disorder is simply a consequence of the rigid band 
theorem 39 : it states that the form of the DOS in an al- 
loy system does not change with alloying, other than via 
a simple translation as given by first-order perturbation 



theory. In our case the magnitude of this shift is given 
by 

AE = (uJ2 c I c p)- xU > ( 33 ) 
v 

where an average over disorder is implied. We can 
confirm that the exact numerical results satisfy quan- 
titatively this expectation by inspection of the data in 
Fig. 14(a). There we plot the position of the minimum 
in the DOS, Ed, for several U and x, being evident that, 
for the concentrations analyzed, the relation (33) is quite 
accurately satisfied up to U ~ 3. For local potentials 
higher than U — 3 the rigid shift of Ed breaks down and, 
in fact, the position of Ed becomes slightly ill defined. 
We point out that concurrently with the shift of Ed (and 
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the band), there is a marked increase in the DOS at Ed, 
unlike the single impurity case (cfr. Fig. 11). This, again, 
is expected and appears already in approximate methods 
like the CPA approximation 30 . Nonetheless, whereas Ed 
shifts linearly for moderate potential strengths, the po- 
sition of the resonance does not vary significantly with 
concentration, and is only enhanced with an increasing 
number of impurities (Fig. 13(b)). 

Another noteworthy aspect of this model has to do 
with the impurity band that emerges at high energies. 
Besides the effects just described, a change in the con- 
centration of impurities implies a concomitant redistribu- 
tion of spectral weight between the main band and the 
impurity band. This is plainly shown in Fig. 14(b) which 
displays the spectral weight in the impurity band against 
the concentration of impurities. This spectral weight is 
calculated by integrating the DOS in the region [4t, oo[. 
As the figure shows, for U > 5t the spectral weight of 
the impurity band saturates at the value x, signaling the 
detachment of the impurity states from the main band. 
For those cases the spectral weight coincides with the 
concentration x. It certainly had to be so because with 
increasing U the impurity band drifts to higher energies, 
eventually disappearing from the problem in the unitary 
limit. As discussed previously in section II A, the spec- 
tral weight of the main band is decreased by precisely x, 
in the presence of a concentration of vacancies of x. This 
is totally consistent with the fact that the local impu- 
rity interpolates between the clean case and the vacancy 
limit. 

Finally, is also clear how the vacancy limit (U — > oo) 
emerges from the data in Fig. 13 as the resonance ap- 
proaches E = and becomes more sharply defined. At 
the same time, the impurity band is displaced toward 
higher and higher energies, eventually projecting out of 
the problem in the vacancy limit. 



D. Non-Diagonal Impurities 

Another effect expected with the inclusion of a sub- 
stitutional impurity in the graphene lattice is the modi- 
fication of the hoppings between the new atom and the 
neighboring carbons. This happens because the host and 
substituting atoms have different radii, because the na- 
ture of the orbitals involved in the conduction band is 
different, or, most likely, a combination of both. Custom- 
ary impurities in carbon allotropes are nitrogen, working 
as a donor, and boron, working as an acceptor 40 . In fact, 
the selective inclusion of nitrogen and/or boron impuri- 
ties in carbon nanotubes is a current practice in the hope 
to tune the nanotubes' electronic response 41 ' 42 ' 43 . 

In general the study of a perturbation in the hopping 
is much less studied in problems with impurities than the 
case of diagonal, on-site, perturbations. In the context 
of our investigations, the perturbation in the hopping 
can, again, be interpreted as an interpolation between 
a vacancy and an impurity. To be more precise, let us 
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FIG. 15: (Color online) Effect of a single substitutional im- 
purity in the LDOS. In panel (a) we plot the LDOS calculated 
at the site of the impurity for the four different values of to 
indicated in each frame. In (b) the situation is identical but 
the LDOS is calculated at the nearest neighboring site of the 
impurity. 



introduce the relevant Hamiltonian: 



H 



c p c P+S 



h.c. . 



(34) 



1,5 



In this case, only nearest neighbor hopping is consid- 
ered. Without the second term, H above is the Hamil- 
tonian for pure graphene. The last sum is restricted to 
the impurity sites, p, and to represents a perturbation 
in the hopping amplitude to its neighbors. Is plain to 
see that, when to = t, all the impurity sites turn into 
vacancies since the hopping thereto is zero. As a result 
of that, this model provides another type of interpolation 
between pure graphene and diluted graphene. An impor- 
tant difference is that this model can be disordered when 
the impurities are placed at random, without breaking 
particle-hole symmetry, and, in this sense, is qualitatively 



14 




Perturbation to the hopping (t Q ) 

FIG. 16: (Color online) Variation of the energy correspond- 
ing to the peak in the LDOS with the magnitude of to- The 
LDOS in question is the LDOS calculated at the impurity 
site. 



much different from the case of local disorder discussed 
in the preceding section. 

We first look at the LDOS in Fig. 15, which contains 
typical results for the local DOS near the impurity, and at 
the impurity site itself. Irrespective of whether the LDOS 
is calculated at or near the impurity, the resulting curves 
display a strong resonance in the low energy region, no 
bound states are formed and the curves are symmetrical 
with respect to the origin. As to increases from zero, 
two simultaneous modifications in these resonances take 
place. The first is that they are clearly enhanced as to 
approaches t. The second is its shift in the direction of 
the Dirac point, in such a way that, when to = 0.9t, the 
peak is already very close to E = 0. With regard to this 
last point, we systematically investigated the variation 
of the peak position in the LDOS at the impurity site 
with the value of to. This dependence, which can be 
seen in Fig. 16, is approximately linear and, for to > 0.6, 
is reasonably well approximated by the linear function 
e max ~ t — to. The apparent saturation for smaller to is 
due to the proximity to the van Hove singularity. The 
study of a single substitutional impurity has been also 
undertaken in ref. 44, with identical results. 

The double-peak structure close to the Dirac point can 
be qualitatively understood from the results regarding a 
vacancy. Suppose that one completely severs the hop- 
ping between a given atomic orbital and its immediate 
neighbors (i.e: set to = t). In this case we are left with 
an isolated orbital with energy E = and a vacancy in 
the honeycomb lattice, which we know also has a zero 
energy mode. If now to is changed slightly, it will cause 
the hybridization of the two zero energy modes with the 
consequent splitting of the energy level, and hence the 
double peaked structure of the LDOS close to the Dirac 
point. 

When we go from one impurity to a finite density of 
impurities, x, we obtain a measurable influence in the 
thermodynamic limit. Our method in this case, consists 
in placing impurities at random positions in the lattice, 




FIG. 17: (Color online) The DOS corresponding to the model 
Hamiltonian (34), with a finite density of impurities. The 
three panels correspond to different values of the perturb- 
ing hopping (to = 0.5, 0.8, 0.9 and 0.95), and within each 
panel the three curves were obtained at different concentra- 
tions (x = 0.01, 0.05 and 0.1). The inset of the bottom panels 
is a magnification of the region near E — 0. 



keeping their concentration constant. The global DOS, 
averaged over several realizations of disorder, is presented 
in Fig. 17. For intermediate values of to, the perturba- 
tion in the hopping induces a resonance appearing at 
roughly the same energies as the ones found in Fig. 16. 
The resonance is enhanced at higher concentrations of 
impurities, and becomes more sharply defined as to — > t. 
Interestingly, as can be seen in the last panels of Fig. 17 
and its inset, the resonant peak splits at higher pertur- 
bations. This splitting depends on the concentration of 
impurities being more pronounced for larger concentra- 
tions and is a new feature introduced by the finite number 
of impurities. As happened already in the case of local 
impurities, the exact numerical results have qualitative 
and quantitative features that could not be anticipated 
from calculations with a single impurity within the usual 
approximation methods. We would also like to point out 
the fact that, from inspection of the above figures, the 
DOS remains zero at E = 0, notwithstanding the sharp 
resonances in its vicinity. Since this model of disorder 
interpolates between clean graphene and graphene with 
vacancies, we are led to a situation similar to the one en- 
countered in sec. II B for uncompensated vacancies. As 
before, it seems that, as the vacancy limit is approached, 
the DOS remains zero at the Fermi energy, despite di- 
verging arbitrarily close to this point, and so the ques- 
tion of the DOS exactly at E = for vacancies lingers. 
Furthermore, unlike what happens with local impurities, 
there is no impurity band nor any high energy features 
appearing as to — > t: the action is all on the low energy re- 
gions. Strictly speaking, in the limit to = t, the impurity 
sites become isolated from the carbon network. Hence 
those sites have to be removed from the Hilbert space for 
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a meaningful physical description of the vacancy case as 
the limit to — > t (for local impurities the removal of the 
impurity sites is akin to the drift of the impurity band 
to infinity, carrying the spectral weight associated with 
the number of the impurities, which projects out of the 
problem). 

Before closing, just a comment on the physical origin 
of this perturbation. In effect, the presence of a substi- 
tutional impurity like N or B will introduce, simultane- 
ously, a perturbation in the hopping, and in the local 
energy. However, it is more or less clear from the discus- 
sions in the previous section that the clearest resonances 
near Ep occur when the local potential, U, is moderate 
or high, which is not the case for boron or nitrogen sub- 
stituents. Hence, the perturbation in the hopping should 
perhaps be more significant in dictating the changes in 
the low-energy electronic structure in the real physical 
system. 

III. CONCLUSIONS 

In this paper we have studied the influence of local 
disorder in the electronic structure of graphene, within 
the tight-binding approximation of eq. (1). We focused 
on vacancies in an otherwise perfect graphene plane and 
the not so extreme cases of local (diagonal) impurities 
and substitutional (non-diagonal, or both) impurities. In 
all cases we saw that disorder brings dramatic alterations 
of the spectrum in the vicinity of the Fermi level. This 
is highly significant since many of the peculiar physical 
properties of graphene stem from the vanishing of the 
DOS at the Dirac point. 

In the case of vacancies, the DOS features a strong di- 
vergence at and close to E = 0, which is associated with 
the formation of quasi- localized states decaying as ~ 1/r 
around the vacancies, which remain even in the presence 
of next-nearest neighbor hopping. Rather interesting is 
the particular case of lattices with uncompensated va- 
cancies, in which case we found the appearance of a gap 
at low energies proportional to the concentration x, and 
the coexistence of localized zero modes in the middle of 
this gap. For the extreme limit of dilution among sites 
of a given sublattice only, we showed that the gap is ro- 
bust, and that a macroscopic number of quasi- localized 
zero modes dominates the spectral density in the middle 
of the gap. Moreover, these zero modes are strictly non- 
dispersive as imposed by symmetry, and give a contribu- 
tion xS(E) to the gapped DOS. This is very interesting, 
in particular if one reasons in terms of magnetic insta- 
bilities and formation of local magnetic moments. Such 



states might be at the origin of local magnetic moments, 
which would explain the magnetism seen experimentally 
in the irradiation experiments 21 . 

We showed how the vacancy case emerges as the lim- 
iting case of a local impurity. In this case the exact cal- 
culation with a single impurity problem was presented, 
taking into account the full dispersion of the honeycomb 
lattice. The results of approximate methods such as CPA 
were subsequently compared with the exact numerical 
solution of the problem with finite concentrations of im- 
purities, and we identified the values of the parameters 
for which these approximations qualitatively break down. 
The discussion of non-diagonal impurities provided yet 
another alternative view of the interpolation between 
clean graphene and vacancies, with relevance for systems 
with dopants that replace the host carbon atoms in the 
honeycomb lattice. One important aspect of the results 
with a finite concentration of these impurities regards the 
splitting of the low energy peaks (insets of Fig; 17), which 
is not captured at a single particle level. The effect has 
to do with situations in which substitutional impurities 
appear close to each other, causing interference and hy- 
bridization effects that lead to the re-splitting of the low 
energy resonances. 

Finally, the results provided for the DOS and LDOS 
are directly testable in real-life samples through scanning 
tunneling spectroscopy techniques and, moreover, the ef- 
fects on the global DOS should reflect themselves in the 
electric transport. For example, one might be able to 
distinguish whether the main effect of a substitutional 
impurity occurs through the modification of the hopping 
to its neighbors, or through the introduction of a local 
potential. 

Note added — The results described here have been 
originally presented in reference 26 during 2006. While 
preparing this manuscript we became aware of the 
preprint 45 with some overlapping results regarding local 
impurities. 
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